library(metafor)
library(plyr)
library(tidyverse)
library(magrittr)
library(foreign)
library(estimatr)
library(emmeans)
library(broom)
library(countrycode)


t4 <- "https://github.com/thomasjwood/kz_multicountry/raw/main/t4.rds" %>% 
  url %>% 
  gzcon %>% 
  readRDS %>% 
  group_by(
    geo_country_code, issue
  ) %>% 
  nest

# estimating models -------------------------------------------------------
t4$mods.w1 <- t4$data %>% 
  map(
    \(i)
    
    broom::tidy(lm_robust(
      out_num ~ treatment,
      data = i
    )
    ))

t4$mods.w2 <- t4$data %>% 
  map(
    \(i)
    
    broom::tidy(lm_robust(
      out_num.w2 ~ treatment,
      data = i
    )
    ))

# overall

w1 <- t4 %>% 
  unnest(mods.w1) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

w2 <- t4 %>% 
  unnest(mods.w2) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

overall.w1 <- rma(w1$estimate, sei = w1$std.error)
overall.w2 <- rma(w2$estimate, sei = w2$std.error)

# global 1

w1 <- t4 %>% 
  filter(issue == "global 1") %>%
  unnest(mods.w1) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

w2 <- t4 %>% 
  filter(issue == "global 1") %>%
  unnest(mods.w2) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

global1.w1 <- rma(w1$estimate, sei = w1$std.error)
global1.w2 <- rma(w2$estimate, sei = w2$std.error)

# global 2

w1 <- t4 %>% 
  filter(issue == "global 2") %>%
  unnest(mods.w1) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

w2 <- t4 %>% 
  filter(issue == "global 2") %>%
  unnest(mods.w2) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

global2.w1 <- rma(w1$estimate, sei = w1$std.error)
global2.w2 <- rma(w2$estimate, sei = w2$std.error)

# country specific

w1 <- t4 %>% 
  filter(issue == "country specific") %>%
  unnest(mods.w1) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

w2 <- t4 %>% 
  filter(issue == "country specific") %>%
  unnest(mods.w2) %>%
  filter(term == "treatmentMisinformation & Fact-check") %>%
  rowid_to_column(., "id")

cs.w1 <- rma(w1$estimate, sei = w1$std.error)
cs.w2 <- rma(w2$estimate, sei = w2$std.error)

# plot

t5 <- bind_rows(tidy(overall.w1), tidy(overall.w2),
                tidy(global1.w1), tidy(global1.w2),
                tidy(global2.w1), tidy(global2.w2),
                tidy(cs.w1), tidy(cs.w2)) %>% bind_cols(meta = 
                                                          c("Overall", "Overall",
                                                            "Gates Foundation  (Global)", 
                                                            "Gates Foundation  (Global)",
                                                            "DNA Modification (Global)", "DNA Modification (Global)",
                                                            "Country Specific", "Country Specific"),
                                                        wave = rep(c("T1", "T2"), 4)) %>%
  mutate(wave = forcats::fct_relevel(wave, "T2"))


t5 %>% 
  mutate(
    meta = meta %>% 
      factor(
        c("Country Specific",
          "Gates Foundation  (Global)",
          "DNA Modification (Global)", 
          "Overall")
      ),
    wave = wave %>% 
      factor(
        str_c("T", 1:2)
      )
  ) %>% 
  ggplot() +
  geom_hline(
    yintercept = 0,
    size = .4,
    linetype = "dashed"
  ) +
  geom_linerange(
    aes(
      ymin = estimate %>% subtract(std.error %>% multiply_by(1.96)),
      ymax = estimate %>% add(std.error %>% multiply_by(1.96)),
      x = wave
    )
  ) +
  geom_point(
    aes(
      wave, estimate
    ),
    shape = 21,
    fill = "white",
    size = 11
  ) +
  geom_text(
    aes(
      wave, estimate, 
      label = estimate %>% 
        round(2) %>% 
        str_replace("0.", ".") %>% 
        str_c(
          c("***", "**", "*", "") %>% 
            extract(
              findInterval(
                p.value, 
                c(-Inf, .001, .01, .05, Inf)
              )
            )
        )
    ),
    family = "Roboto",
    size = 2.5
  ) +
  facet_wrap(
    ~meta, nrow = 1
  ) +
  scale_x_discrete(
    labels = 1:2
  ) +
  labs(x = "Wave",
       y = "Meta-Analysis Estimate")